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Abstract. 

The extended effective multiorbital Bose-Hubbard-type Hamiltonian which takes 
^ ■ into account higher Bloch bands, is discussed for boson systems in optical lattices. 

, It is shown that the renormalization of Hamiltonian parameters depends on the 

dimension of the problem studied. Therefore, mean field phase diagrams do not scale 
with the coordination number of the lattice. The effect of Hamiltonian parameters 
renormalization on the dynamics in reduced one-dimensional optical lattice potential 
■ is analyzed. We study both the quasi-adiabatic quench through the superfluid-Mott 

CSJ . insulator transition and the absorption spectroscopy, that is energy absorption rate 

when the lattice depth is periodically modulated. 
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1. Introduction 

Ultracold bosonic atom gases in optical lattices have been an ultrahot research area 
in recent years (for the recent reviews and an extensive reference list see [IH3|). They 
provide means to create and control experimental systems mimicking different condensed 
matter physics models [3H6]. The interest has been stimulated in part by the fact that 
there exists [7j an accurate mapping of continuous Hamiltonian to a Hamiltonian on 
a lattice — the Bose- Hubbard (BH) Hamiltonian, originally formulated by Gersch and 
Knollman [5J. The lattice models significantly ease the analytical [9TTT2] and numerical 
analysis, though promising new ideas were recently proposed that can enable the analysis 
in continuous variables |13] beyond the mean field level [TJJ[T5]. 

Bosons in one dimensional (ID) optical lattices have also been the area of 
extensive experimental research [T6T[2~1~] . The corresponding ID BH Hamiltonian can be 
effectively addressed numerically by Density Matrix Renormalization Group (DMRG) 
and related techniques [22]. These techniques have broad applications, in particular 
enable simulations targeting real-life many-body systems [23], with controllable error 
and no systematic errors due to unsounded assumptions. Their success relies on area 
laws that control the growth of entanglement [24J. The entanglement is used as a 'small 
parameter' [25 ] |26 ] that makes it possible the construction of an efficient variational 
set - the so called Matrix Product States (MPS). Several numerical investigations of 
experimental and 'close to experimental' systems have been performed [271 - 129] . They 
focused largely on two aspects: a quench through a phase transition and the simulation 
of absorption spectroscopy [16]. The phase transition from the superfluid to the Mott 
insulator phase occurs when the lattice depth is increased beyond a critical value [U|30]. 
The adiabaticity of this process has been addressed in [T5ll29p 31 | . In the second example, 
the energy absorption rate is analyzed as a function of the frequency of modulation of 
the lattice depth [271132, 33J. The locations, number of peaks in the spectrum and their 
heights are related to the state of the gas: either superfluid or insulating. 

The derivation of the BH Hamiltonian assumes the restriction of physics to the 
lowest Bloch band of the optical lattice. This assumption is reasonable, as the energy 
gap between the first and second Bloch bands is in most cases around 10 times larger 
than the energy scale in the discrete model. It is also much larger than the thermal 
energy scale k^T providing additional argument for a zero temperature analysis. An 
early analysis |3l] suggested that higher bands may be included by an appropriate 
modification of BH Hamiltonian parameters for large occupation numbers. The effects 
due to higher Bloch bands have been also studied for lower densities [35T438] in relation 
to the discovery that higher bands affect the superfluid-insulator transition in Bose- 
Fermi mixtures [391B0] . This gave explanation for experimental observations of the shift 
of the Superfluid (SF) - Mott insulator (MI) phase transition [41] which could not be 
explained by single-band approaches. In [36] , the first excited band is included in a two- 
flavour model; effects of higher bands could also be built in via an effective three-body 
interaction in the lowest band [35|l3~7]. 
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Recently, another approach for 3D optical lattices has been proposed |3l? | BS | H3" ] 
(which somehow resemble in spirit |3l] while being used for moderate atomic densities). 
In this approach, the higher bands are included in the on-site Hamiltonian, which is then 
diagonalized and many-body ground states of this problem for different total number 
of particles are used as a local Hilbert space, replacing the usual Fock basis \n). The 
multiband lattice Hamiltonian is then expressed in this new basis, yielding an effective 
single-band model with occupation dependent parameters — renormalized values of the 
initial Hamiltonian parameters for the lowest Bloch band. Interestingly, the change of 
BH Hamiltonian parameters is large in the MI regime, where the energy gap is also 
large, contrary to a naive intuition. 

In this paper, we consider the derivation of the BH model's effective coupling 
constants, taking a closer look at the underlying, quite challenging, numerical problem 
- the accurate diagonalization of the onsite Hamiltonian. We study how the 
dimensionality of the optical lattice affects the renormalization scheme. The dependence 
on dimension implies that the mean field diagrams of the system no longer depend solely 
on the coordination number of the lattice. We then take a look at the effects of the 
renormalization of coupling constants on the dynamics in a ID optical lattice: both 
quench through the SF-MI transition and the absorption spectroscopy are analyzed. We 
find that renormalization of atom-atom interactions significantly shifts and sometimes 
modifies the absorption peaks. It also appears to have a serious effect on adiabaticity 
predictions for the SF-MI phase transition. 



2. Tight-binding descriptions of an ultracold boson gas in an optical lattice 

The second quantization Hamiltonian for a dilute gas of bosonic atoms interacting by 
s-wave scattering (with scattering length a s ) in the optical lattice potential V(r) and 
external trapping potential V e (r) is of the form: 



H 




+ V e (r)jijj(r) + — ^(f) , 0^(r)'0(r)'0(r) — /i?/^(r)?/>(r) 



dr{\) 



with g = 4irh 2 a s /m. The Hamiltonian admits a natural energy scale, set by the 
recoil energy Er = where k — ^ and A denotes the optical lattice wavelength. 

We assume this energy unit from now on. The optical lattice potential is typically: 
V(x, y, z)/Er = s x cos 2 {kx) + s y cos 2 (ky) + s z cos 2 (kz). If s := s x <C s y = s z =: s±, then 
a classical setup for a ID optical lattice is obtained |16) . Tunneling in y, z directions 
is highly suppressed and the system may be considered as a series of independent ID 
tubes along the x direction. Similarly when s := s x = s y <C s z =: s± a 2D optical 
lattice is obtained. When s x = s y = s z =: s the potential corresponds to a 3D lattice. 
Let us emphasize that all calculations of renormalized parameters presented below are 
truly three-dimensional ones, the labels "ID" or "2D" just refering to situations where s 
and s± are very different. 
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Jaksch and Zoller in their seminal paper [7J introduced a mapping of (OQ) onto the 
BH Hamiltonian (here Vi is the local energy shift due to V e : Vi — V e (fi)): 

H B h = - oaa) + h.c. + — £ h^fii - 1) - £ - K). (2) 

The field operator is expanded in the set of the lowest band Wannier functions of the 
lattice: ip(x) = Y^ w i{ x ) a i, ~ ~ here wf denotes the (real valued) Wannier function 
localized on site i of the lattice for the a'th Bloch band (we shall denote sites by roman 
subscripts and bands by greek superscripts). The parameters U and J are expressed by 
the appropriate integrals of Wannier functions: J = — J w®{f) ( — + V{rfj Wj(r) 
(where i and j are neighbouring sites) and U = 47r ^ as / Wi(f) 4 dr. 

Without restriction to the lowest Bloch band, the expansion in the full Wannier 
functions basis set would give a multiband variant of the Bose-Hubbard Hamiltonian: 

a,i,j a...S 

i...l 

+ J2(E a + V i -fi)(a^at (3) 

a.i 

with 
and 

By construction J°f = for a ^ (3. For sufficiently deep optical lattices (typically 
s > 2), it is enough to restrict hopping to nearest-neighbor sites, as tunneling amplitudes 
are exponentially damped with the hopping distance (for shallow lattice next nearest 
neighbours hopping may be necessary - see |21]). 

A 3D Wannier function being a product of ID Wannier functions, the 3D integral 
in U^] 6 is a product of 3 integrals over each coordinate. The interaction parameters 
differ for each direction, as Wannier functions depend on the lattice depth, which may 
be different in each direction. 

The Hamiltonian ([3]), although mathematically exact, is difficult to use in practice, 
even in ID systems, because the onsite dimension d of the n-particle problem restricted 
to the lowest B bands is d = ~ t n n_1 ), (the onsite problem is genuinely 3D even for 
a quasi-lD models) and moreover the numerical complexity of the best ID algorithms 
scales with d at least as 0(d 3 ). The complexity of more sophisticated approaches (such 
as MERA, PEPS [4H447] ) is several orders higher. 

Thus for computational purposes, one must restrict the local Hilbert space. 
Assuming that interactions are on-site only, i.e. U^] 6 ^ for i = j = k = I together 
with considering the lowest Bloch band only (a = /3 = 7 = <5 = 0) [7] leads directly to 
the Bose-Hubbard Hamiltonian, Eq. (j2J), provided we chose the zero of the energy axis 
at Eq. 
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A more sophisticated approach is discussed in |39[I43|. The on-site Hamiltonian, 
restriction of to a single site: 

Hioc = H E + H U = J2 E a h a + U^aiala^as. (4) 

a 0/875 

(with n a = aJ a a a ) can be diagonalized to yield a space of n particle ground states. 
The eigenenergies eg of the on-site n particle ground states IV'o)) are the starting point 
in determining new values of U parameters in the effective Hamiltonian. To define 
renormalized values of U, the energy has to be decomposed into the interaction 
energy [which in case of the BH model, is just ^n(n — 1)] and a single-particle energy 
(which in the BH case shifts \i by the lowest Bloch band energy). The most natural way 
to define the interaction energy would be to use: 

^n(n - 1) = (CI E U^fafaatW). (5) 

a/375 

Unfortunately, U n cannot be defined in such a way if we request the Hamiltonian to 
have a form resembling Eq. (j2J). That is because the single particle energy is no longer 
a linear function of n. This can be circumvented by defining U n via: 

= Y n ( n + nE °- (6) 

This definition makes Hamiltonians (j2j) and ((Tj) similar. But U n is no longer the 
interaction energy, it also contains contributions of higher Bloch bands population to 
the single-particle energy. From now on we use definition (jBJ). Note that U n depends 
nontrivially on the geometry of the lattice. 

The second stage is to reintroduce inter-site couplings. Even if only one band 
is taken into account (so that U n is simply U), the inter-site interaction 
induces an effective coupling, which is proportional to the density, that is a term 
Uitij° a i a j( n i + n j ~ 1) + ^- c - m ^ ne effective Hamiltonian (called bond-charge term 
in [39 | I48|). Similarly to the original tunneling term, it is important only between 
nearest neighbors as soon as s is larger than unity. In the low-density regime (typically 
n < 7), this contribution leads to an increase (since < for s > 1) of the tunneling 
amplitude J — > J — U^j (ni + rij — 1). The correction is at most of the order of the 
raw tunneling. When higher bands are taken into account, the modification of |?/>q) 
induces a renormalization of the standard tunneling term (as well as of the bond-charge 
term) which becomes also dependent on the occupation numbers of sites between which 
tunneling occurs. 

The effective multiorbital (EMO) Hamiltonian finally becomes: 

h emo = _ £ + hmCm) + £ E^ n{n _ i) P i, ( 7 ) 

n.i,rij n,i 

where P % n = \i,n)(i,n\. 

This Hamiltonian will allow us to study the influence of higher bands on the 
dynamics later. First we discuss the accurate numerical determination of the U and 
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J parameters that, in itself, is a challenging problem, giving an insight into the physics 
involved. 

2.1. Solving the onsite problem 

The single site problem is equivalent to finding the n particle ground state of the 
Hamiltonian (Jl]). If the lowest B Bloch bands are included (typical values of B : 4 
[42 j, 9 - |39]). the problem quickly becomes too involved computationally to be 
exactly diagonalized, and truncation of the basis has to be performed. 

For small particle numbers, the problem is dominated by the He term. The 
direct transition from the lowest to the first excited band is forbidden by symmetry 
consideration so the relevant energy scale is A = E<i — Eq. Promoting a particle from 
the lowest to the second excited band, via interaction term is proportional to n 3 / 2 lj 2000 
where [J 2000 is the corresponding interaction integral. This yields for the validity 
of the perturbative approach the condition n 3 / 2 lj 2000 <g; A. For typical parameters 
corresponding to Rb scattering length and s of the order of 30, A/U 2000 ~ 60 yielding 
the limiting value of n <C 15. 

For small occupation numbers, a perturbative approach seems justified. 
Perturbation theory enables us to estimate the impact each excited vector \ip p ) has on 
the ground state energy. Let |^o) be an n-boson ground (Fock) state of He- The larger 
the matrix element \(ipo\Hioc\' l l , p)\ 2 and the smaller the energy (ippl'Hiod'ipp) the larger is 
the impact. The perturbative scheme provides a hint on how to choose an "optimal" 
subset of the basis in which the full problem could be diagonalized. Since the exact 
diagonalization in the variational basis is the last step, we do not follow perturbation 
theory exactly, but only its spirit, just to choose a close to optimal basis, not to calculate 
the energy correction. Details of basis generation and variational space sampling are 
given in the Appendix. A more traditional method [39] is to choose a subset according 
to least energy principle — with minimal (Vvl^ZodVv) ■ 

We have performed a detailed analysis comparing both methods for s = s± = 34.8, 
a strongly coupled case (3D optical lattice, 2a s /X = 0.014, A = 754 nm - parameters 
taken as typical values from [39]). We choose a system with n = 2 — 5 particles and 
40000 basis vectors according to both least-energy (as in [39J) and the perturbative 
method (39900 vectors are generated within the first order, and 100 within the second 
order perturbative scheme). 

The ground state energy - obtained from numerical diagonalization of the on-site 
Hamiltonian in a restricted subset - versus the number B of bands included is shown 
in Figure [H for various n. The least energy method clearly gives the false impression 
of saturation of results if B « 7 — 9 bands are included. The false saturation occurs 
because the least-energy method does not evaluate \(ipo\Hioc\4 , p)\ 2 an d therefore fills 
the variational basis with low-energy irrelevant vectors with vanishingly small matrix 
elements. The perturbative-like approach does not show similar saturation effects and, 
moreover, suggests that linear extrapolation of the results may be performed. We find 
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that the best compromise between computational effort and accuracy is to perform 
extrapolation of results as a function of 1/B by means of the ansatz U n (B) = U^° + Cq/B. 
The same extrapolation scheme can be used for the J parameters (but leads to less 
drastic modification of the results'). 
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Figure 1. Comparison of the effective on-site interaction strength U n obtained using 
diagonalization of the on-site Hamiltonian on two different basis sets with the same 
size equal to 40000. If basis vectors are chosen according to their energy (dashed lines), 
false saturation effects appear. Estimating the influence by a perturbative-like scheme 
(solid lines) does not seem to suffer from saturation effects. The 3D case is considered: 
s = sj_ = 34.8. 



The U n parameters for the ID, 2D and 3D lattices are presented in Figure [21 while 
the renormalized tunneling amplitudes are shown in Figure [3j We find that, in low 
dimensions, U n vary less with s compared to the full 3D lattice. The high transverse 
lattice causes significant renormalization of U n even for small s, as s± is still large. 

Inspection of Fig. [1] and Figj2] shows that the difference between consecutive U n is 
approximately constant, i.e., 

tf„-i -U n *W, (8) 

at least for low densities (for typical lattice parameters, W is constant up to 10%). This 
is easily understood: the alternative effective theory of [35] expresses the correction to 
the on-site energy term via a three-body interaction term (Eq. (12) of [35]) and W is 
simply related to their parameter U3. The deviations from the linear form, eq. (jSJ), 
may be then related to higher order terms in [35], i.e. four-body term, etc. Similarly, 
the lowest order perturbative term discussed above gives a correction to the interaction 
energy term n{n — l)U/2 of the order n 3 [77 2000 ] 2 /A. That yields a crude estimate for 
W w U 2000 /A w 1/60, in good agreement with Fig. CD 
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Lattice depth s 

Figure 2. Renormalized interaction parameters U n vs. strength of the optical lattice, 
for various dimensionalities (black solid — ID, red dashed — 2D, blue dot dashed — 
3D). Interaction and lattice parameters: 2a s /X = 0.014, A = 754 nm. The transverse 
lattice height is sj_ = 34.8. The curves meet at s = s±. 




5 10 15 20 25 5 10 15 20 25 30 
Lattice depth s 



Figure 3. Renormalized tunneling amplitudes J ni .n 2 vs - strength of the optical lattice, 
for various dimensionalities (black solid — ID, red dashed — 2D, blue dot dashed — 
3D). Interaction and lattice parameters: 2a s /\ = 0.014, A = 754 nm, s± = 34.8. The 
curves meet at s = sj_. 
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3. Consequences of coupling constants renormalization 

3.1. Mean field diagrams for different lattice dimensions 

A homogeneous (without external trap, Vi = 0) Hamiltonian may be taken to the 
thermodynamic limit. Then the particle density is determined by the chemical potential 
fi and the Hamiltonian parameters: J and U. A MI phase is determined by integer 
density (n) and noncompressibility: = 0. The rest of the phase diagram is the 

compressible SF phase [HI EJ HHJ ED]- We perform now the mean field analysis of the 
phase diagrams of Hamiltonians fl2J) and (jTJ). A Gutzwiller analysis of the Bose- Hubbard 
model is a variational minimization of the following functional (i.e. mean ground state 
energy) : 

H BH [^\ = (i;\H BH m = -ZzJiaiK&j) + j(Hn ~ 1)) - (9) 

using the Gutzwiller ansatz: = (g) with \ipi) = fn\ n ) the on-site wavefunction. 
The influence of the lattice geometry is reduced only to the coordinate number, z, in 
the first term, as (a*) = (aj) due to the translational invariance in thermodynamic limit. 
The phase diagram depends only on the single parameter 

For the EMO Hamiltonian, the dependence of interaction parameters on the 
dimensionality of the optical lattice is nontrivial. We shall use the data from Figure [2j 
Thanks to the translational invariance, the Gutzwiller mean field approach to the 
effective Hamiltonian (J3) is equivalent to the minimization of the following functional: 



H emo [4j] = -2z £ J ni ,nM\a>i)(n 1 ^i)mal\n 2 )(n 2 \4j l ) + 

ni,n2 

+ E - l)l<«| 2 - ^£™l<^>| 2 - (10) 

n n 

Clearly a single parameter is no longer sufficient to describe the mean field problem. Let 
us denote by Jbh the tunneling amplitude and by Ubh the interaction of the standard 
BH hamiltonian. As is strictly decreasing with the lattice height s, it provides 

a way to plot results calculated for lattice depth s in the ( "ubh ' ^ coordinate space 
used for a traditional phase diagram. This mapping allows also to directly compare 
the results obtained using Eqs.(jSJ) and ffTU]) . Figure H] shows the Gutzwiller phase 
diagrams for ID, 2D, 3D lattices. In contrast with the ordinary BH model, there is 
a nontrivial dependence on the dimension. Let us stress again that phase diagrams for 
the MO parameters J nunj /U n do not need to be directly related to Z V ^ H . The physical 
observable that is common for the ordinary Bose-Hubbard phase diagram and the MO 
effective theory is the lattice depth s. Moreover, for dimensions 1 and 2, there is actually 
a whole family of different phase diagrams indexed by s±. We show just a single choice 
for a generic value of s± = 34.8. 

In ID, the mean field approximation is inaccurate. To get a reliable phase diagram, 
we have used energy minimization through imaginary time evolution using the Time 
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Evolving block Decimation (TEBD) [25|[2l3] algorithm. We fix the lattice size to be 
L = 100 (we have checked that choosing a larger lattice size L = 200, 300, 400 does 
not alter the results significantly, except at MI tips, where an approximate finite size 
scaling is performed). This is significantly less computationally demanding than using 
the infinite, translationally invariant version of TEBD [471[51]- The transverse lattice 
height is again s± = 34.8. Let us denote by E(N, s) the ground state energy of 
a iV-particle system for lattice height s. We calculate approximations to the critical 
values of chemical potentials delimiting a Mott insulator region with average filling n by 
~ E(nL + 1, s) — E(nL, s), /i_(s) ~ E(nL, s) — E(nL — 1, s). We plot the phase 
diagrams for both the BH and the EMO Hamiltonians in Figure [5j We again see similar 
results: the Mott lobes shrink also in ID, as predicted by the mean field approach. As 
shown below, this is also reflected in the dynamical properties. 

3 
2.5 

2 

as 

§ 1.5 

1 

0.5 



"0 0.05 0.1 0.15 0.2 

zJ BH (s)/U BH (s) 

Figure 4. Mean field phase diagrams for ID, 2D, and 3D lattices. Different curves 
denote borders between MI and SF phases. Dashed black lines correspond to the 
standard BH model for any dimension, red, green, blue curves denote ID, 2D and 
3D lattices of the EMO Hamiltonian. The limit zJbh{s)/Ubh{s), s — > oo is different 
for each dimension. The s — > oo limit corresponds to the ill-defined situation in 
which the transverse lattice is shallower than the main lattice (this formal limit is 
also dimension-dependent). The perpendicular lattice depth is fixed at s± = 34.8, 
A = 754 nm, 2a s /A = 0.014 as appropriate for 87 Rb [39]. 

3.2. Modulation of optical lattice - absorption spectroscopy 

By periodically modulating the lattice depth, one transfers energy to the atomic 
sample in the lattice. Absorption spectroscopy - also incorrectly nicknamed modulation 
spectroscopy - consists in studying the dependance of the energy absorption rate with the 
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Figure 5. ID phase diagram obtained using imaginary time evolution and the TEBD 
algorithm. Black dashed curves presents the standard BH ID case, the red solid lines 
are obtained for the EMO ID model with s±_ = 34.8. 



modulation frequency. This absorption is sensitive to the quantum phases present in the 
system, as shown in early experiments j 16 j . It has been simulated |27" t 1321133] for atoms 
in an optical lattice in the presence of a harmonic confinement, using a standard ID BH 
model. It seems interesting to see whether excited bands affect the absorption spectra. 
To this end, we consider the real time evolution of the ground state of a given system 
at s — So, exposed to a time- varying lattice height s(t) = Sq + s m cos ut, s m /so <C 1. 



The simulation is performed in the presence of an harmonic trap Vi 



K[l 



! . The 



energy of the system is measured after some fixed time. The energy gain (per particle) 
as a function of the modulation frequency yields the absorption spectrum. 

For a deep optical lattice in a predominantly Mott insulating phase, the absorption 
spectrum for the standard BH Hamiltonian consists of a few peaks located at 
multiplicities of U [23 E2J [33] . The situation is slightly more complex for the EMO 
Hamiltonian. The position of peaks can be easily determined in the deep Mott regime 
(J — > 0). States excited during the modulation are mainly those that differ from the 
ground state by nearest-neighbours transfer of one particle [32|I52]. Up to a small 
correction due to the difference of the local chemical potential /ij = fi — Vi, the excitation 
energy is determined by the occupation numbers of the source site i and destination site 
j. It is: 

AE(m] Tij) = - - 

(1 



-U ni ni(ni 



-U nj rij(rij - 1) + U nj+1 nj(nj + 1) 
i) + Un^im - i)(m - 2) + A/Xjj] 



with Afiij = ^ — fij = Vj — Vi. Nearest neighbours excitation means that 



3\ 




Figure 6. Effects of higher Bloch bands on absorption spectroscopy in the deep 
Mott (low J) regime, s = 15, s± = 40. Panel (a) shows the well-known wedding 
cake structure with n = l,n = 2,n = 3 Mott plateaus. Excitations within each 
plateau (colored respectively light gray, dark gray, black, for n = 1,2,3) have energies 
depending on the Mott plateau density and the trapping potential. Inward and 
Outward hopping lead to a splitting of the absorption structure, a partial splitting 
for moderate harmonic trap [(b), k = 0.001] or a broad well resolved structure for 
a shallow trap [(c), k = 0.0001] in contrast to the standard BH case (d). 



and \rii — rij\ < 1. If U n — U (BH Hamiltonian), we have that AE{rii\nj) = 
(nj—rii+l)U +A/ijj. By virtue of Eq. ([8]) we may approximate Eq. (fTTj) by: AE(nf, rij) = 
(nj — ni + l)U ni + (ni—nj — 2)(n.i+nj — l)W. For a trapped gas with maximum occupation 
number n = 3, the relevant values are: AE(1; 1) = U 2 , AE(1; 2) = 3U 3 -U 2 pa 2U 2 -3W, 
AE(2; 2) = 3U 3 - 2U 2 pa U 2 - 3W, AE(2; 3) = 6t/ 4 - 3U 3 - U 2 pa 2U 2 - 9W, 
AE(3;3) = QU 4 - QU 3 + U 2 pa U 2 - 6W : AE(n;n - 1) = 0, with A/i neglected for 
clarity 

A qualitative comparison of the expected absorption spectra for standard BH 
case and the EMO model is possible. The density profile in the large s, low 
hopping, limit shows the well known wedding cake structure (see Figure |6ti). Weak, 
periodic modulation leads mainly to nearest neighbour excitations between any pair or 
neighbouring sites. For the standard BH Hamiltonian, the excitation spectrum consists 
of a large peak at energy U (Figure [6H). The nonzero width of the peak is due to 
variations of the local chemical potential (the presence of a trap): the shallower the 
trap, the narrower the peak. 

For the EMO Hamiltonian, the particle-hole excitations from different Mott 
plateaus have different mean excitation energies. The shift with respect to the mean 
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value is determined by the A/Zjj. For shallow traps A/ijj m 2(j — i)n{i — i ) and the 
excitation spectrum from Mott plateau n consists in two bands — one corresponding to 
Inward (I) hopping, A/ijj < 0, the other one to Outward (O) hopping, A/ijj > 0, with 
respect to the trap centre. For the central Mott plateau with density nmax a single, 
broad peak in the excitation spectrum emerges. This is clearly visible in Figures [6b, c. 
Two cases have been studied: a system of N = 260 particles in a trap with k = 0.001 
[(b), moderate case] and N = 700, k = 0.0001 [(c), very shallow trap]. Both cases were 
studied for the ID optical lattice with s = 15, s± = 40, A = 830 nm, a = 5.1 nm. The 
tunneling J n ^ nj is artificially set to (deep Mott regime). 
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Figure 7. Absorption spectrum (modulation time t = 100Ti/Er, modulation 
amplitude s m = 1). Black dashed lines correspond to the BH model, red solid curves 
to the effective multiorbital theory. Left panel shows spectra on a common energy 
scale, observe the significant shift of the EMO structure toward smaller energies. Bars 
above the plot give the mean expected positions of peaks for the n=2 Mott plateau. 
Right panel shows the same data with rescaled energy axes (Ubh for the black curve, 
U2 for the red one). 

A similar, but smaller system is analysed in a subsequent numerical study, for the 
same parameters taking fully into account the tunneling effects. We choose a much 
tighter trap with curvature, k = 0.009, and use true values of the hopping constants 
Jn t ,nj- The Wannier function calculations give Ubh — 0.662i?#, and the renormalization 
procedure gives U 2 ~ 0.565-Er, W ~ 0.0125-Er. 

We fill the trap with N = 36 particles. This system is similar to the one studied 
in [32|l33], The system states are represented by MPS vectors and evolved using the 
TEBD algorithm |25ll26] . The ground state of the system, being the initial state for the 
evolution is calculated using an imaginary time evolution with bond dimension x — 50. 
The local Hilbert space assumes maximal filling of 6 bosons per site. 

The density profiles of ground states of the BH and EMO models are practically 
the same with a central plateau of 2 particles per site. Thus, any change of excitation 
frequencies can be interpreted as an effect of coupling constants renormalization. We 
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have performed the absorption spectroscopy simulation for time t = lOOk/E^. The 
modulation amplitude of the lattice was s m = 1. The results are presented in Figure [71 

The major difference between the spectra obtained for the BH and EMO 
Hamiltonians is a significant shift of the observed structures. While for the BH case, the 
main structure is centered at Ubh, it has a similar shape, but centered around U 2 in the 
EMO case. Because of the steep harmonic trap - thus the large changes in local chemical 
potentials - the structure of the peaks is rather complex. Note the global broadening for 
the EMO case, and an additional peak in the main U2 structure, corresponding to the 
n — 1 plateau excitations having an energy larger by roughly 3W, as discussed above. 

A second small peak on the right appears at E — 1.75U 2 (EMO case) and 
E = 1.8Ubh (BH case). It corresponds to a particle- hole excitation on the edge between 
the n = 1 and n = 2 Mott plateaus as identified in [33]. The right panel shows that the 
spectra becomes quite similar if rescaled by their proper energy scale, Ubh or U 2 . 

The absorption spectra are quite sensitive to the details of the system. Taking 
the same parameters for a slightly larger number of particles may create situations 
where the ground states of the BH and EMO Hamiltonians significantly differ. This is 
then reflected in the absorption spectra. If the density profile contains a n = 3 or higher 
plateau - then the structure of peaks may become more complicated, as discussed above. 

We have also compared the absorption spectra in the superfluid regime. The lattice 
height is fixed at s = 5, s± = 40 [18J. The system is modulated for t = 50H/Er with 
s m = 0.2. The results are presented in Figure [HI Unlike in the Mott regime, the positions 
of the absorption peaks are not determined solely by the interaction. In particular, no 
global shift of the structure is observed. In both cases, one observes a broad resonance 
around the recoil energy, with complicated detailed structures. Note that the modulation 
depth is much smaller than in the Mott insulator situation (to avoid significant excitation 
of the system) and therefore the absorbed energy per particle is much smaller than in 
Fig. [7J 

3.3. Florence experiment revisited 

In the Florence experiment [18], the starting point is a ultra-cold gas in a harmonic 
trap (without optical lattice). The optical lattices are then ramped up (assuming 
s{t) = 0.2s±(t)) with an exponential ramp s(t) ~ s (l — exp(t/r)), for r = 30 ms, 
and a total ramping duration 100 ms. The system soon becomes quasi-lD producing 
a set of ID tubes. If a "disordered" system is desired, an additional optical lattice, 
with different wavelength A2, is superimposed along the tubes. This adds a potential 
Vd{%) = s 2 sm 2 (k 2 x). For s 2 <C s, it acts effectively as the shift of the on-site energy, i.e. 
an additional pseudo-random disorder AVi = s 2 sin 2 (-^fti + <fj , where <p represents the 
offset between the two optical lattices. We take a generic value <fi = 0.12345. If the ratio 
A/A2 is chosen irrational enough, the bichromatic lattice simulates a disorder well enough 
for a finite system [53ti57] . The dependence of the pseudo- disorder amplitude s 2 on time 
is set by demanding that s 2 (t) ~ s(t) (all optical lattices are ramped up simultaneously). 
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Figure 8. Absorption spectra in the superfluid case. The absorbed energy per particle 
is plotted as a function of the frequency of modulation of the lattice depth. Here 
s = 5, s m = 0.2, sj_ = 40, Ubh ~ 0.465-Er. The renormalized interaction parameter: 
U2 ~ 0A06Er. The left panel corresponds to 20 particles, the right panel to 36 particles. 

We will consider 3 cases: no disorder — 0), weak disorder (s 2 = 55s), and strong 

disorder (s 2 = J^)- 

Consider first the no disorder case S2 = 0. As initial state, we take 151 particles on 
81 sites in the presence of a harmonic confinement coming both from the trap and the 
transverse laser profile. The detailed procedure using the TEBD algorithm is described 
in [291133]. 

After the optical lattice is ramped, absorption spectroscopy is performed for 30ms 
(the conversion unit is 20. 91ft/ Er = 1ms). In the recent numerical investigation of 
this experiment [29], a discrepancy between experimental [H] and numerical results was 
found. The reported position of the first absorption peak was 1.9 kHz [18j, while Wannier 
function calculations gave 2.3kHz [29j. The renormalization procedure renormalizes the 
value of the U parameter to U 2 — 2 kHz, U3 = 1.96 kHz, U4 = 1.91 kHz. This 
suggests that the positions of absorption peaks due to the n = 1, 2, 3 Mott plateaus are: 
U 2 = 2kHz, -2U 2 + 3U 3 = 1.85kHz, U 2 - 6U 3 + 6U4 = 1.74kHz. The "average" peak 
position is 1.87 kHz. Therefore, the EMO Hamiltonian provides an estimate of the peak 
position in good agreement with the experiment. 

A simulation of absorption spectrum performed for this system confirms this finding 
as shown in Figure [9j Although the initial state when the periodic lattice modulation 
starts is not the ground state, but a wavepacket dynamically created during the ramping 
of the lattice, the peak positions are well predicted by the EMO model. The position of 
the first and second peaks agree quite well with the experiment (the relative height is 
different presumably because of the strong modulation used in [18]) . 

The exponential ramping of the optical lattice in the experiment [18] may not be 
adiabatic as discussed in |29| using a standard BH description. Due to the discrepancy 
in the position of the absorption peak, the lattice depth was adjusted in |29j. Instead of 
ramping the lattice up to s = 16, the final value s — 14 was considered. In some sense, 
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Figure 9. Absorption spectrum obtained by applying lattice modulation with 
amplitude s m = 1 on the wavepacket created by exponential ramp up to s = 16. 
The black dashed line corresponds to the standard BH model, the red line is the result 
of effective multiorbital theory. The position of the absorption peaks in the latter case 
reproduce well the experimental results |18| . 



such a simplified approach may be viewed as a renormalization of the BH parameters 
(without insight into its origin explained in section l3~2l) . Let us stress that the agreement 
between the experimental position of absorption peaks and the EMO predictions prove 
the necessity of using the effective multiorbital theory to explain quantitatively the 
experimental results. 

With that modification of the final s value, it was found using the BH model |29] 
that the overlap of the prepared wavepacket on the ground state at the final s value 
was about 9% in the absence of disorder. It is most interesting to see how taking into 
account higher bands within effective multiorbital theory affects the adiabaticity of the 
dynamics. The simulation performed for the similar exponential ramp starting at s = 5 
up to s = 16 yields an overlap of the dynamical wavepacket on the ground state at 
s — 16 equal to 17.3%. It may be qualitatively understood: in the EMO model, the 
effective interactions are weaker and the effective tunneling larger allowing particles to 
redistribute more efficiently among sites during the ramp. 

We have also tested an optimized s(t) pulse shape as in [29] - By choosing s(t) 
changing slowly close to the phase transition point for the n = 3 Mott lobe, we have 
been able to enhance adiabaticity up to 33% squared overlap with the ground state. 

The presence of disorder has devastating influence on adiabaticity, similarly to the 
standard BH case [29]. We have found that for a small disorder, S2 = -^s the squared 
overlap is a fraction of percent (0.005) while for the strong disorder S2 = -^s it becomes 
vanishingly small (of the order of 10~ 9 , beyond the accuracy of the calculation). 
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4. Conclusion 

The aim of this paper is two-fold. In the first part, we have presented an efficient 
numerical implementation of the approach sketched in [39] which makes it possible to 
compute the parameters of the effective Hamiltonian for bosons in optical lattices. The 
approach goes beyond the standard Bose-Hubbard model [7] limited to the lowest Bloch 
band. The effective Hamiltonian approach which includes contributions from higher 
lying bands (multiorbital approach) has been shown to lead to new effects and even new 
phases |37] for bosonic systems (see also |35j). Our scheme of perturbatively generated 
basis seems clearly superior to the energy-selected basis used in [39] for low and moderate 
occupation numbers and allows for better estimates of Hamiltonian parameters. These 
estimates may be extrapolated to an infinite number of bands. 

We have applied the method not only to the standard 3D cubic lattice, but also to 
reduced ID and 2D problems, where the lattice depth is different in various directions. 
The effective Hamiltonian obtained depends on the dimensionality of the problem. In 
effect, mean field phase diagrams as obtained with the Gutzwiller ansatz, differ even 
if they are rescaled by the lattice coordination number. It turns out that the role of 
excited bands is even more pronounced for reduced dimensionality problems than for 
a 3D lattice. 

Motivated by this difference, we have investigated whether the dynamics is different 
in a standard Bose-Hubbard model and for the effective multiorbital theory. We have 
considered two cases, the energy absorption created by modulation of the lattice height 
and the quasi-adiabatic passage from the superfluid to the Mott insulator phase. In 
both situations, it turns out that taking into account the density dependent tunneling 
terms as well as modification of interactions may lead to significant differences between 
two approaches. For the same lattice depth, the effective interactions turn out to 
be significantly weaker than in the standard Bose-Hubbard model. This results in 
profound differences in the absorption spectra such as significant shifts of absorption 
peaks. Similarly, the full effective theory predicts that the transition from superfluid to 
Mott insulator is more adiabatic than with the standard Bose-Hubbard model |29j. 

The results presented in the present paper should have a direct applicability to any 
experiment using bosons in an optical lattice, with multiple site occupations. 
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Appendix A. Diagonalization of the onside Hamiltonian 

We describe the approach we use to generate a perturbatively based variational set used 
in the diagonalization of the onsite n particle problem: 

U loc = H E + Hu = J2 E a n a + £ U a ^ s alamos. (A.l) 

The most elementary possible excitation promotes a boson from the a = to the 
a = 2 band. This defines two limits: n 3 / 2 U 2000 < E 2 - E and n 3 / 2 U 2000 » E 2 - E . 
Ijap-yS Decome smaller as a, f3, 7, S indices grow. 

The first regime corresponds to the ordinary Bose-Hubbard model in which Hu can 
be treated as a small perturbation, with zero order ground state = (aj ^ |0), just 
as assumed for the ordinary BH hamiltonian. The large density regime is dominated by 
Hjj. A straightforward application of the multiband model in that regime is not justified, 
however a mean field based renormalization scheme leads to the approximate model of 
the same form with modified parameters |34j . The transition to the nonperturbative 
regime occurs at approximately 10-15 atoms per site. 

0.9 1 ' 1 1 1 1 1 1 1 1 1 1 1 1 1 ' 1 1 I 




Number of particles n 



Figure Al. Comparison of U n obtained from diagonalization of T-Li OCl Eq. dA.ll) . 
with energy selected basis and with perturbatively chosen sets. Each perturbative 
set consists of 10000 vectors of first order and of 0, 100, 800 vectors of second order 
(curves: black, red, green). Blue dashed curve denotes the least-energy basis result. 
For small and moderate n, the perturbative based basis is clearly superior over the 
least-energy set. The failure of perturbation theory approach for n sufficiently large is 
apparent too. 



Consider the low energy regime. By vectors reachable in k-ih order perturbative 



expansion, we call those Fock states 
basis can be generated with order 



tp p ) for which (i/; \H^\'t/)p) 7^ 0. In particular a full 
. Let us denote by the set of vectors reachable 
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in k-th order and unreachable in k — 1-th order. The full variational basis of Fock states 



is B = U £>fc, with Bk pairwise disjoint, and B^ = 0, for k > 

k 

For numerical diagonalization, a proper, not too large S C B has to be chosen. 
In [39], S consists of vectors \ip) with the least = (ip[Hioc\' l P) • I n the first order 
perturbation method, we choose vectors from B\ with the largest values of: 



/lW = ln JM^m. (A . 2) 

Disregarding vectors from U B& seems to be a crude approach, but still our numerical 

k>2 

calculations prove this "basic" perturbative approach to be more efficient for low 
densities than the least-energy-based selection. The perturbation theory provides 
a way to evaluate a perturbative contribution of vectors from Bk, for arbitrary k, 
which, unfortunately, is computationally involved (summation over intermediate states, 
degeneracy resolution). As we do not calculate the ground state energy within the 
perturbation theory treatment, but only motivate the choice of a variational basis for 
numerical diagonalization, the following function: 



/ 2 (V0=hi sup — — — — r. (A.3) 

|Vi}eSi y^\4>) ~ ^\4>o)){^\ipi) ~ ^W'o)) 

is chosen to approximate the relevancy measure for vector It minimizes the state's 
\xfj) energy and maximizes the overlap. 

Detailed comparison of the three methods: the least-energy, first order perturbative 
and "improved" first order perturbative approach is presented in Figure [TJ We have 
used K = 10000 vectors (plus additional second higher order vectors for the "improved" 
method), set s = s± = 38 and compared the three methods as a function of n. Two 
regimes: perturbative and not perturbative emerge, as expected. In the small and 
moderate n regime, the perturbative approach gives better estimate for the ground 
state energy than the least-energy method [39]. The opposite tendency is visible in the 
nonperturbative regime. 

Let us describe in detail how a choice of the basis with the largest f\ and / 2 values 
is performed. We have fixed the maximal number of Bloch bands included at B = 15. 
We use Markov chain Monte Carlo method which is quite general and can be applied at 
any order k of perturbation expansion. B^ is in general too large to evaluate function 
fk for all elements (its size increases exponentially with the total number of particles). 
It is usually possible just to scan the whole B\ set, so from this point on we assume 
k > 2. To choose K vectors with largest f^ values, we construct a random walk based 
on Metropolis' Monte Carlo algorithm. A state of the random walk is a finite k + 1-tuple 
of n-particle Fock states V = (l^o); l^i); ^2)1 ■ ■ ■ 1 l^fci); 1"^)); f° r l^i) £ ^t- For any such 
V we define the following generalization of (IA.3I) : 

™ , \{MHu\^i){MHu\^)\. ■ ■ {^ k -i\Hu 
9k{y) — tn 



(E w - E^ o) )(E\^ k _ 1 ) - E\^ o) ) . . . {E\^) - E\^ o) ) 
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sup 

i/>i>eB 1) ...,|Vfe-i>6Bfc_i 



9k(V). 



(A.4) 



To get the random walk, we have to update V. First we choose at random 
a Fock state 6 V to be updated. With equal probability, we update one or two 
particles of preserving the total parity of the state. One particle update is done 



normalized. Direction y is not special in any way: with equal probability any of x, y, z 
is chosen. After the update a proposition V' is prepared. We automatically reject 
updates for which |^) ^ B\. If that is not the case, the acceptance probability is 
determined as in Metropolis algorithm: it is given by min{l, exp[/3(g%(V) — fi'fc(V))]}. 
The inverse temperature f3 is tuned to optimize sampling efficiency — we choose it by 
requiring the acceptance rate to be close to 0.3. After a successful update, the last 
element of the tuple V, state \ipk) is accepted into the solution set if its perturbative 
importance <?fc(V) is in the K lowest values recorded so far. The accepted vector \ipk) 
is memorized as well as the importance value g^iV 1 ). If \ipk) had been generated before, 
the memorized value of g^ is updated (only if the new value of larger than the old one). 
If, in a subsequent few thousand sweeps (empirical value), no vector makes it into the 
solution set, nor gj, values are updated, then the procedure is restarted. This ensures 
that all low energy excitations are taken into account (the starting point is always the 
low energy configuration). Altogether, we make 2 x 10 9 MC sweeps to generate basis of 
size 40000 (as used for the results presented in the main part of the paper). 

If all Bloch bands were included, then the set of Fock space would be infinite. On 
the other hand, only a finite number of them could satisfy the inequality: f^ > s. The 
values of fi for the remaining states are very close to 0, and a singularity in density 
of states arises. Logarithm is used to "smoothen" this singularity for numerical 
purposes. It does not affect the ordering, as In is increasing injection. 
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